Bond-order modulated staggered flux phase for the t— J model on the square lattice 
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Motivated by the observation of inhomogeneous patterns in some high-T c cuprate compounds, sev- 
eral variational Gutzwiller-projected wave-functions with built-in charge and bond order parameters 
are proposed for the extended t — J — V model on the square lattice at low doping. First, following a 
recent Gutzwiller-projected mean-field approach by one of us (Phys. Rev. B. 72, 060508(R) (2005)), 
we investigate, as a function of doping and Coulomb repulsion, the stability of the staggered flux 
phase with respect to small spontaneous modulations of squared unit cells ranging from 2 x 2 to 
\/32 x \/32- It is found that a 4 x 4 bond-order (BO) modulation appears spontaneously on top of the 
staggered flux pattern for hole doping around 1/8. A related wave-function is then constructed and 
optimized accurately and its properties studied extensively using an approximation-free variational 
Monte Carlo scheme. Finally, the competition of the BO-modulated staggered flux wave-function 
w.r.t. the d-wave RVB wave-function or the commensurate flux state is investigated. It is found 
that a short range Coulomb repulsion penalizes the d-wave superconductor and that a moderate 
Coulomb repulsion brings them very close in energy. Our results are discussed in connection to the 
STM observations in the under-doped regime of some cuprates. 

PACS numbers: 74.72.-h,71.10.Fd,74.25.Dw 



I. INTRODUCTION: MODELS AND METHODS 

The observation of a d-wave superconducting gap 
in the high-T c cuprate superconductors suggests 1 that 
strong correlations are responsible for their unconven- 
tional properties and superconducting behavior. The 
two-dimensional (2D) t— J model is one of the simplest 
effective models proposed 2 to describe the low energy 
physics of these materials, 

Ht-J = ~* S ('•-'•- ' 

+J]TS 4 -S, (1) 

(id) 

The electrons are hopping between nearest neighbor sites 
of a square lattice leading to a kinetic energy term (first 
term of^) as well as an exchange energy due to their spin 
interaction (second term), where Si denotes the spin at 
site i: Si = \c\ a a a ^Ci^ and a is the vector of Pauli ma- 
trices, stands for a pair of nearest neighbors. H t -j 
operates only in the subspace where there are no dou- 
bly occupied sites, which can be formally implemented 
by a Gutzwiller projector (see later). In the following we 
set \t\ — 1 (unless specified otherwise) and we adopt a 
generic value of t/J = 3 throughout the paper. Because 
of the particle-hole symmetry in the square lattice the 
sign of t does not play any role. Although this model is 
formulated in a very simple form, the nature of the quan- 
tum correlations makes its physics very rich, and even the 
ground state of the t—J hamiltonian was not yet charac- 
terized for finite doping and large cluster size. However, 



the t—J model was investigated extensively by unbiased 
numerical techniques 3 - as well as by mean-field* and vari- 
ational Monte-Carlo approaches^. All approaches found 
a d-wave superconducting phase and a phase diagram 
which accounts for most of the experimental features of 
the high-T c cuprates*^. In the limit of vanishing doping 
(half- filling) , such a state can be viewed as an (insulat- 
ing) resonating valence bond (RVB) or spin liquid state. 
In fact, such a state can also be written (after a simple 
gauge transformation) as a staggered flux state (SFP)*^, 
i.e. can be mapped to a problem of free fermions hopping 
on a square lattice thread by a staggered magnetic field. 

Upon finite doping, although such a degeneracy 
breaks down, the SFP remains a competitive (non- 
superconducting) candidate with respect to the d-wave 
RVB superconductor. In fact, it was proposed by 
P.A. Lee and collaborators^ 1 ** 2 ** 3 - that such a state bears 
many of the unconventional properties of the pseudo- 
gap normal phase of the cuprate superconductors. This 
simple mapping connecting a free fermion problem on 
a square lattice under magnetic field** to a correlated 
wave-function (see later for details) also enabled to con- 
struct more exotic flux states (named as commensurate 
flux states) where the fictitious flux could be uniform and 
commensurate with the particle densitji****. In this par- 
ticular case, the unit-cell of the tight-binding problem is 
directly related to the rational value of the commensurate 
flux. 

With an increasing number of materials and novel 
experimental techniques of constantly improving reso- 
lution, novel features in the global phase diagram of 
high-T c cuprate superconductors have emerged. One of 
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the most striking is the observation, in some systems, 
of a form of local electronic ordering, especially around 
1/8 hole doping. Indeed, recent scanning tunnelling mi- 
croscopy/spectroscopy (STM/STS) experiments of un- 
derdoped I^S^CaC^Os+a (BSCO) in the pseudogap 
state have shown evidence of energy-independent real- 
space modulations of the low-energy density of states 
^T3Qg^i7 i i8 i i9 i 2o w it n a spatial period close to four lat- 
tice spacings. A similar spatial variation of the electronic 
states has also been observed in the pseudogap phase of 
Ca 2 _ 2; Na :r Cu02Cl2 single crystals (x = 0.08 ~ 0.12) by 
similar STM/STS techniques^. Although it is not clear 
yet whether such phenomena are either generic features 
or really happening in the bulk of the system, they nev- 
ertheless raise important theoretical questions about the 
stability of such structures in the framework of our mi- 
croscopic strongly correlated models. 

In this paper, we analyze the stability and the proper- 
ties of new inhomogeneous phases (which may compete in 
certain conditions with the d-wave superconducting RVB 
state) by extending the previous mean-field and varia- 
tional treatments of the RVB theory. In addition, we 
shall also consider an extension of the simple t— J model, 
the t— J— V model, containing a Coulomb repulsion term 
written as, 

V=lj2 V (\ i -i\)(m-n)(n j -n) , (2) 

where n is the electron density (N e /N , N e electrons on 
a iV-site cluster). Generically, we assume a screened 
Coulomb potential : 



where we will consider two typical values £q = 2, 4 and 
Vo € [0, 5] and where the distance r is defined (to mini- 
mize finite size effects) as the periodized distance on the 
torusS^. The influence of this extra repulsive term in the 
competition between the d-wave RVB state and some in- 
homogeneous phases is quite subtle and will be discussed 
in the following. 

To illustrate our future strategy, let us recall in more 
details the simple basis of the RVB theory. It is based 
on a mean-field hamiltonian which is of BCS type, 

Hbcs = ^2 (~Xo4r c > + h - c ) 

+ Yl ( A ^ c h c ]i + h - c ) n ^ - ( 4 ) 

i,<7 

where xo is a constant variational parameter and Ajj is 
a nearest neighbor d-wave pairing (with opposite signs 
on the vertical and horizontal bonds) and /i is the chem- 
ical potential. As a matter of fact, the BCS mean-field 
hamiltonian can be obtained after a mean-field decou- 
pling of the t—J model, where the decoupled exchange 



energy leads to the xo and Ajj order parameters. In 
this respect, we expect that the BCS wave- function is a 
good starting point to approximate the ground state of 
the t—J model. However, such a wave-function obviously 
does not fulfill the constraint of no-doubly occupied site 
(as in the t—J model). This can be easily achieved, at 
least at the formal level, by applying the full Gutzwiller 
operator— Vg = lli(l — n »T n u) to tne BCS wave-function 
I^BCs): 

IV'RVb) = Vg IV'bcs) • (5) 

The main difficulty to deal with projected wave- functions 
is to treat correctly the Gutzwiller projection Vg. This 
can be done numerically using a conceptually exact vari- 
ational Monte Carlo (VMC) techniqu o 5 i 6 i 7 on large clus- 
ters. It has been shown that the magnetic energy of 
the variational RVB state at half-filling is very close 
to the best exact estimate for the Heisenberg model. 
Such a scheme also provides, at finite doping, a semi- 
quantitative understanding of the phase diagram of the 
cuprate superconductors and of their experimental prop- 
erties. Novel results using a VMC technique associated to 
inhomogeneous wave- functions will be presented in Sec- 
tion EH 

Another route to deal with the Gutzwiller projection 
is to use a renormalized mean-field (MF) theory^ in 
which the kinetic and superexchange energies are renor- 
malized by different doping-dependent factors gt and gj 
respectively. Further mean-field treatments of the inter- 
action term can then be accomplished in the particle- 
particle (superconducting) channel. Crucial, now well 
established, experimental observations such as the exis- 
tence of a pseudo-gap and nodal quasi-particles and the 
large renormalization of the Drude weight are remarkably 
well explained by this early MF RVB theory 8 . An exten- 
sion of this approaches will be followed in Section E] 
to investigate inhomogeneous structures with checker- 
board patterns involving a decoupling in the particle-hole 
channel. As (re-) emphasized recently by Anderson and 
coworkers^, this general procedure, via the effective MF 
Hamiltonian, leads to a Slater determinant |^mf) from 
which a correlated wave- function Vg \^mf) can be con- 
structed and investigated by VMC. Since the MF ap- 
proach offers a reliable guide to construct translational 
symmetry-breaking projected variational wave- functions, 
we will present first the MF approach in scctionlTTlbefore 
the more involved VMC calculations in Section ITTT1 

II. GUTZWILLER-PROJECTED MEAN-FIELD 
THEORY 

A. Gutzwiller approximation and mean-field 
equations 

We start first with the simplest approach where the 
action of the Gutzwiller projector Vg is approximately 
taken care of using a Gutzwiller approximation scheme^. 
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We generalize the MF approach of Ref. l25l to allow for 
non-uniform site and bond densities. Recently, such a 
procedure was followed in Ref. to determine under 
which conditions a 4 x 4 superstructure might be stable 
for hole doping close to 1/8. We extend this investigation 
to arbitrary small doping and other kinds of supercells. 
In particular, we shall also consider 45-degree tilted su- 
percells such as \/2 x y/2, y/E x y/8 and \/32 x V32- 

The weakly doped antiferromagnet is described here 
by the renormalized t— J model Hamiltonian, 



mz = -*<* E (4<^> + h - c -) + j 9j E s * 



(6) 



where the local constraints of no doubly occupied sites 
are replaced by statistical Gutzwiller weights gt — 
2x/(l + x) and gj = 4/(1 + x) 2 , where x is the hole 
doping. A typical value of tj J = 3 is assumed hereafter. 

Decoupling in both particle-hole and (singlet) particle- 
particle channels can be considered simultaneously lead- 
ing to the following MF hamiltonian, 

HmF = -t E 9ij(cl,a C 3,v + h - C -) + E eini ' a 
(ij)tr io 

\> E slMiA° c w + " \Xij\ 2 ) ( 7 ) 

li e <^A,' : :,<i , • '>■<■■ x, h. 

where the previous Gutzwiller weights have been ex- 
pressed in terms of local fugacities z$ = 2xi/(l 
Xi) (xi is the local hole density 1 — (rii)), 



and g 



J _ 



ah 



(2 — Zi)(2 — Zj), to allow for small 



non-uniform charge modulations 2 -. The Bogoliubov-de 
Gennes self-consistency conditions are implemented as 

Xji = ( cj j,a c i : cr) &nd Aj'j = {Cj,—<rCi,<r) — {ci, — cr c j,cr)- 

In principle, this MF treatment allows a description 
of modulated phases with coexisting superconducting or- 
der, namely supcrsolid phases. Previous investigations 26 
failed to stabilize such phases in the case of the pure 2D 
square lattice (i.e. defect-free). Moreover, in this Sec- 



tion, we will restrict ourselves to Aj 



0. The case 



where both and \ij are non-zero is left for a future 
work, where the effect of a defect, such as for instance a 
vortex, will be studied. 

In the case of finite Vq , the on-site terms ej may vary 
spatially as — /i+ej, where fj, is the chemical potential and 
ej are on-site energies which are self-consistently given by, 



E^ 



(8) 



In that case, a constant J2i^j + " 2 ) nas to 

be added to the MF energy. Note that we assume here a 
fixed chemical potential \x. In a recent work 2 ^, additional 
degrees of freedom where assumed (for Vq — 0) imple- 
menting an unconstrained minimization with respect to 
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FIG. 1: 4 x 4 unit-cell used in both the MF approach and 
the variational wave-function. Note the existence of 6 inde- 
pendent bonds (bold bonds) , that for convenience are labelled 
from 1 — 6, and of 3 a priori non- equivalent sites. The center 
of the dashed plaquette is the center of the (assumed) dv 
symmetry. Other sizes of the same type of structure are con- 
sidered in the MF case, respectively : 2 x 2, v8 x v8, and 
\/32 x \/32 unit cells. 



the on-site fugacities. However, we believe that the en- 
ergy gain is too small to be really conclusive (certainly 
below the accuracy one can expect from such a simple 
MF approach). We argue that we can safely neglect the 
spatial variation of /j, in first approximation, and this will 
be confirmed by the more accurate VMC calculations in 
Section lnTl Incidently, Ref. |22| emphasizes a deep connec- 
tion between the stability of checkerboard structures 2 ^ 
and the instability of the SFP due to nesting properties 
of some parts of its Fermi surface 2 ^. 



B. Mean- field phase diagrams 

In principle, the mean-field equations could be solved 
in real space on large clusters allowing for arbitrary 
modulations of the self-consistent parameters. In prac- 
tice, such a procedure is not feasible since the number 
of degrees of freedom involved is too large. We there- 
fore follow a different strategy. First, we assume fixed 
(square shaped) supercells and a given symmetry within 
the super-cell (typically invariance under 90-degrees rota- 
tions) to reduce substantially the number of parameters 
to optimize. Incidently, the assumed periodicity allows us 
to conveniently rewrite the meanfield equations in Fourier 
space using a reduced Brillouin zone with a very small 
mesh. In this way, we can converge to either an absolute 
or a local minimum. Therefore, in a second step, the MF 
energies of the various solutions are compared in order 
to draw an overall phase diagram. 

In previous MF calculations 26 , stability of an inhomo- 
geneous solution with the 4x4 unit-cell shown in Fig. 
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was found around x — 1/8. Here, we investigate its sta- 
bility for arbitrary doping and extend the calculation to 
another possible competing solution with a twice-larger 
(square) unit-cell containing 32 sites. The general so- 
lutions with different phases and/or amplitudes on the 
independent links will be refered to as bond-order (BO) 
phases—:. Motivated by experimentsi£>21, a C^v sym- 
metry of the inhomogeneous patterns around a central 
plaquette will again be assumed for both cases. Note 
that such a feature greatly reduces the number of varia- 
tional parameters and hence speeds up the convergence 
of the MF equations. Starting from a central plaquette, 
like in Fig. ^ a larger V32 x V32 unit-cell (not shown) 
can easily be constructed with 10 non-equivalent bonds 
(with both independent real and imaginary parts) and a 
priori 6 non-equivalent sites. Note that this new unit-cell 
is now tilted by 45 degrees. 

At this point, it is important to realize that patterns 
with a smaller number of non-equivalent bonds or sites 
are in fact subsets of the more general modulated struc- 
tures described above. For example, the SFP is obvi- 
ously a special case of such patterns, where all the Xi,j 
are equal in magnitude with a phase oriented to form 
staggered currents, and where all the sites are equiva- 
lent. This example clearly indicates that the actual struc- 
ture obtained after full convergence of the MF equations 
could have higher symmetry than the one postulated in 
the initial configuration which assumes a random choice 
for all independent parameters. In particular, the equi- 
librium unit-cell could be smaller than the original one 
and contain a fraction (1/2 or 1/4) of it. This fact is 
illustrated in Fig. [21 showing two phase diagrams pro- 
duced by using different initial conditions, namely 4x4 
(top) and a/32 x a/32 (bottom) unit-cells. Both diagrams 
show consistently the emergence of the SFP at dopings 
around 6% and a plaquette phase (2x2 unit-cell with 
two types of bonds) at very small doping^i*^. In ad- 
dition, a phase with a a/8 x a/8 super-cell is obtained 
for a specific range of doping and Vq (see Fig. [21 on the 
top). Interestingly enough, all these BO phases can be 
seen as bond-modulated SFP with 2, 4, 8 and up to 16 
non-equivalent (staggered) plaquettes of slightly differ- 
ent amplitudes. This would be consistent with the SFP 
instability scenario^ which suggests that the wavevec- 
tor of the modulation should vary continuously with the 
doping. Although this picture might hold when Vq = 0, 
our results show that the system prefers some peculiar 
spatial periodicities (like the ones investigated here) that 
definitely take place at moderate Vq. 

Let us now compare the two phase diagrams. We find 
that, except in some doping regions, the various solutions 
obtained with the 4x4 unit-cell are recovered starting 
from a twice larger unit-cell. Note that, due to the larger 
number of parameters, the minimization procedure start- 
ing from a larger unit-cell explores a larger phase space 
and it is expected to be more likely to converge to the ab- 
solute minimum. This is particularly clear (although not 
always realized) at large doping x = 0.14, where we ex- 



1 - 







2x2 



SFP 



x 
|oc 



4x4 



FL 



0.05 0.1 0.15 0.2 

doping 



1.5 - 



0.5 



2x2 



I N 



SFP 



4x4 



FL 



0.05 



0.1 0.15 

doping 



0.2 



FIG. 2: Mean-field phase diagrams obtained by solving sell- 
consistently the mean-field equations on a 128 x 128 lattice 
(for lo — 4) vs hole doping x and repulsion Vo (in units of 
J). Top: results obtained assuming a 4 x 4 unit cell; Bottom: 
same with a a/32 x a/32 tilted unit cell. In both cases, a d v 
symmetry is assumed (see text). 



pect an homogeneous Fermi Liquid (FL) phase (all bonds 
are real and equal) , as indeed seen in Fig. on the bot- 
tom. On the contrary, Fig. [21 on the top reveals, for 
Vq/ J G [1.5, 3], a modulated a/8 x a/8 structure, which is 
an artefact due to the presence of a local minimum (see 
next). 

Since the MF procedure could accidentally give rise to 
local minima, it is of interest to compare the MF en- 
ergies obtained by starting with random values of all 
independent parameters within the two previously dis- 
cussed unit-cells. For convenience, we have substracted 
from all data either the FL (in Fig.^a)) or the SFP (in 
Fig- [Hlb)) reference energy. From Figs.[3{a,b) we see that 
we can converge towards a local energy minimum, often 
modulated in space, which is not the absolute minimum. 
Indeed, over a large doping range, the lowest energy of 
the all solutions we have found is obtained for homoge- 
neous densities and bond magnitudes. Nevertheless, we 
see that the 4x4 modulated phase is (i) locally sta- 
ble and (ii) is very close in energy to the homogeneous 
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(SFP) phase which, often, has a slightly lower energy. 
Note that, around x ~ 1/8, the states with y/8 x y/8 
and \/32 x v32 supercells are clearly metastable solu- 
tions (and using a larger initial unit-cell is not favorable 
in the latter case). In contrast, in this range of doping, 
the 4x4 checkerboard state is very competitive w.r.t. 
the SFP. Therefore, it makes it a strong candidate to be 
realized either in the true ground state of the model, or 
present as very low excited stated. In fact, considering 
such small energy differences, it is clear that an accurate 
comparison is beyond the accuracy of the MF approach. 
We therefore move to the approximation-free way of im- 
plementing the Gutzwiller projection with the VMC tech- 
nique, that allows a detailed comparison between these 
variational homogeneous and inhomogcncous states. 



III. VARIATIONAL MONTE CARLO 
SIMULATIONS OF 4 x 4 SUPERSTRUCTURES 
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FIG. 3: (Color online) (a) Energy per site (in units of J and 
for t — 3 J) obtained by solving the mean-field equations using 
the initial 4x4 unit-cell (see text) for a moderate value of 
Vo- The SFP energy is also shown for comparison. The FL 
energy has been substracted from all data for clarity, (b) 
Comparison of the energies (for Vo = 0) using different initial 
conditions (see text), a 4 x 4 or a ^32 x V32 unit-cell; due 
to very small energy differences, the SFP energy is used as 
a reference for an easier comparison. The different phases 
specified by arrows and characterized by the number of sites 
Nsc of their actual supercells refer to the ones in Fig.H For 
doping x = 0.14, the minimization leads to a solution with 
small imaginary parts (of order 10 -4 ) very similar to a FL 
phase, which we call FL*. 



Motivated by the previous mean-field results we have 
carried out extensive Variational Monte Carlo simula- 
tions. In this approach, the action of the Gutzwiller pro- 
jection operator is taken care of exactly, although one has 
to deal with finite clusters. In order to get rid of discon- 
tinuities in the d-wave RVB wave- function, we consider 
(anti-)periodic boundary conditions along e y (e x ). As a 
matter of fact, it is also found that the energy is lower for 
twisted boundary conditions, hence confirming the rele- 
vance of this choice of boundaries. We have considered a 
16 x 16 square cluster of N = 256 sites. We also focus on 
the 1/8 doping case which corresponds here to N e = 224 
electrons on the 256 site cluster. Following the previous 
MF approach, we consider the same generic mean-field 
hamiltonian, 

H M F = (- kj c L c 3<y + h - c ) + £ i C ia C iv > ( 9 ) 

where the complex bond amplitudes tij can be written as 
lijj | e l8i - j , and &ij is a phase oriented on the bond 
The on-site terms allow to control the magnitude of the 
charge density wave. However, the energy was found to 
be minimized for all the equal to the same value in the 
range Vq = [0, 5] and for the two parameters lo = 2, 4. In 
fact, we find that strong charge ordered wave- functions 
are not stabilized in this models. 

In this Section, we shall restrict ourselves to the 4x4 
unit-cell where all independent variational parameters 
are to be determined from an energy minimization. This 
is motivated both by experimentsiiSi and by the previ- 
ous MF results showing the particular stability of such a 
structure (see also Ref . |26|) . As mentioned in the previous 
Section, we also impose that the phases and amplitudes 
respect the C±v symmetry within the unit-cell (with re- 
spect to the center of the middle plaquette, see Fig. QJ, 
reducing the numbers of independent links to 6. To avoid 
spurious degeneracies of the MF wave-functions related 
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to multiple choices of the filling of the discrete k-vectors 
in the Brillouin Zone (at the Fermi surface), we add very 
small random phases and amplitudes (10 -6 ) on all the 
links in the 4x4 unit cell. 

Let us note that commensurate flux phase (CFP) are 
also candidate for this special 1/8 doping. In a previ- 
ous study, a subtle choice of the phases (correspond- 
ing to a gauge choice in the corresponding Hofstadter 
problem^) was proposed^, which allows to write the 
4> = p/16 (p < 16) flux per plaquette wave-function 
within the same proposed unit-cell^ and is also expected 
to lead to a better kinetic energy than the Landau gauge 
(in the Landau gauge the unit-cell would be a line of 
16 sites). However, we have found that the CFP wave- 
functions turned out not to be competitive for our set 
of parameters Vo, due to their quite poor kinetic energy, 
although they have very good Coulomb and exchange 
energies. We argue that such CFP wave-functions would 
become relevant in the large Coulomb and/or J regimes 
(see table HJ. 

In order to further improve the energy, we also add a 
nearest-neighbor spin-independent Jastrow^ term to the 
wave- function, 

Vj = exp I a ^ n.nj , (10) 

V / 

where a is an additional variational parameter. Finally, 
since the t — J model allows at most one fermion per site, 
we discard all configurations with doubly occupied sites 
by applying the complete Gutzwiller projector Vg. The 
wave- function we use as an input to our variational study 
is therefore given by, 

IWax) = VgVj IV^Mf) (11) 

In the following, we shall introduce simple notations for 
denoting the various variational wave- functions, BO for 
the bond-order wave function, SFP for the staggered flux 
phase, RVB for the d-wave RVB superconducting phase, 
FS for the simple projected Fermi-Sea, and we will use 
the notation MF/J (MF = BO, SFP, RVB, FS) when 
the Jastrow factor is applied on the mean-field wave- 
function. Finally, it is also convenient to compare the 
energy of the different wave-functions with respect to 
the energy of the simple projected Fermi-Sea (i.e. the 
correlated wave-function corresponding to the previous 
FL mean- field phase), therefore we define a condensation 
energy as e c = e var - e FS . 

In Fig. 0] we present the energies of the three wave- 
functions BO/J, SFP/J and RVB/J for Coulomb po- 
tential Vo E [0,5]. We find that for both £ = 2 and 
^o = 4 the RVB phase is not the best variational wave- 
function when the Coulomb repulsion is strong. The 
bond-order wave-function has a lower energy for Vq > 2 
and £ = 2 (V > 1.5 and £ = 4). Note that the 
(short range) Coulomb repulsion in the cuprates is ex- 
pected to be comparable to the Hubbard U, and therefore 



TABLE I: Set of energies per lattice site for Vo = 1 and £o = 
4 for different wave-functions. The best commensurate flux 
phase in the Landau gauge with flux per plaquette p/16 was 
found for p — 7. We also show the energy of the CFP with 
flux 7/16 written with another choice of gauge. We show the 
total energy per site (E to t), the kinetic energy per site (-Et), 
the exchange energy per site (-Ej) and the Coulomb energy 
per site (Ev). 



wave-function 


Etot 


Et 


E, 


E v 


FS 


-0.4486(1) 


-0.3193(1) 


-0.1149(1) 


-0.0144(1) 


CFP 7/16 1 


-0.3500(1) 


-0.1856(1) 


-0.1429(1) 


-0.0216(1) 


CFP 7/16 2 


-0.4007(1) 


-0.2369(1) 


-0.1430(1) 


-0.0208(1) 


SFP 


-0.4581(1) 


-0.3106(1) 


-0.1320(1) 


-0.0155(1) 


BO 


-0.4490(1) 


-0.3047(1) 


-0.1302(1) 


-0.0141(1) 


RVB 


-0.4564(1) 


-0.3080(1) 


-0.1439(1) 


-0.0043(1) 


SFP/J 


-0.4601(1) 


-0.3116(1) 


-0.1315(1) 


-0.0169(1) 


BO/J 


-0.4608(1) 


-0.3096(1) 


-0.1334(1) 


-0.0177(1) 


RVB/J 


-0.4644(1) 


-0.3107(1) 


-0.1440(1) 


-0.0086(1) 



Landau gauge 
2 Gauge of Ref[3 



TABLE II: Order parameters for the different wave-functions 
for Vo = 1.5 and £o = 4. We depict the following order param- 
eters: tij x e^»'L where tij (4>i,j) is the amplitude (phase) 
of {cfcj}, and the exchange energy (Si.Sj), for the 6 inde- 
pendent bonds labelled for convenience according to Fig. Q 
The sign of 4>i,i is according to the staggered flux pattern 
(see arrows in Fig. |UJ. We note that the RVB/J is uniform 
by construction. The variational superconducting order pa- 
rameter is Aiivb = 0.3 for the RVB/J wave-function and 
Arvb — for the SFP/J and BO/J wave-functions. 





bond 1 


bond 2 


bond 3 


bond 4 


bond 5 


bond 6 


RVB/J 
SFP/J 
BO/J 


0.077(1) 
0.085(1) 
0.082(1) 


0.077(1) 
0.085(1) 
0.083(1) 


0.077(1) 
0.085(1) 
0.093(1) 


0.077(1) 
0.085(1) 
0.088(1) 


0.077(1) 
0.085(1) 
0.086(1) 


0.077(1) 
0.085(1) 
0.084(1) 


\<t>i,j\ 
RVB/J 
SFP/J 
BO/J 




0.438(1) 
0.527(1) 




0.438(1) 
0.502(1) 




0.438(1) 
0.473(1) 




0.438(1) 
0.390(1) 




0.438(1) 
0.338(1) 




0.438(1) 
0.384(1) 


-{Si.Sj) 
RVB/J 
SFP/J 
BO/J 


0.215(1) 
0.197(1) 
0.215(1) 


0.215(1) 
0.197(1) 
0.207(1) 


0.215(1) 
0.197(1) 
0.215(1) 


0.215(1) 
0.197(1) 
0.187(1) 


0.215(1) 
0.197(1) 
0.186(1) 


0.215(1) 
0.197(1) 
0.170(1) 



Vo ~ 5 or 10 seems realistic. Independently of the rela- 
tive stability of both wave-functions, the superconduct- 
ing d-wave wave-function itself is strongly destabilized 
by the Coulomb repulsion as indicated by the decrease of 
the variational gap parameter for increasing Vq and the 
suppression of superconductivity at Vo ~ 7 (see Fig.|S}. 
Nevertheless, we observe that the difference in energy 
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FIG. 4: Energy per lattice site of the RVB/J, SFP/J 
and BO/J wave-functions minus the energy of the projected 
Fermi- Sea wave- function. 




FIG. 5: Kinetic and exchange energy per site of the RVB/J 
wave-function minus the respective exchange and kinetic en- 
ergy of the simple projected Fermi-Sea. Inset: value of the 
variational d-wave gap. 



FIG. 6: Total energy per site of the BO/J minus the energy 
of the SFP/J wave- functions. 




0.1 0.2 0.1 0.2 0.3 0.4 0.5 

<t> 

FIG. 7: Total energy per site of the BO/J variational wave- 
function with variational parameters Im (Uj) = i<p on the 
bonds 1, 2, 3, and Im (t iyj ) = ±0.149 on the bonds 4, 5, 6. The 
sign of Im (tij) is oriented according to the staggered flux 
pattern. We have chosen for all the links Re ( tij) = 0.988. 
Results for Vo = and Vo — 5 with £o — 4 are shown. 



between the bond-order wave-function and the staggered 
flux phase remains very small. We show in table [H] the 
order parameters measured after the projection for the 
RVB/J, SFP/J and BO/J wave-functions. As ex- 
pected the RVB / J and the SFP/ J wave- functions are 
homogenous within the unit-cell. In contrast, the BO/J 
wave- function shows significant modulations (expected to 
be measurable experimentally) of the various bond vari- 
ables w.r.t their values in the homogeneous SFP. In Fig.|H| 
we show the small energy difference (see scale) between 
the two wave- functions. Interestingly, the difference is 
increasing with the strength of the potential. We notice 
that the two wave-functions correspond to two nearby 



local minima of the energy functional at zero Coulomb 
potential (see Fig.[7J), which are very close in energy (the 
BO/J wave- function is slightly lower in energy than the 
SFP/ J) and are separated by a quite small energy bar- 
rier. Note that in Fig. [7] we consider the variational bond 
order parameters and not the projected quantities. 

When the repulsion is switched on, the height of the 
energy barrier increases and the SFP/J wave-function 
does not correspond anymore to the second local min- 
ima. Indeed, when Vo > the second local energy min- 
ima shifts continuously from the point corresponding to 
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FIG. 8: Kinetic, exchange and Coulomb energy per site of the 
BO I J wave-function minus the respective associated energy 
of the SFP/ J wave- function. 



the simple SFPj J wave- function. The metastable wave- 
function lying at this second local minima is a weak bond- 
order (SFP-like) wave-function that preserves better the 
large kinetic energy while still being able to optimize 
better the Coulomb energy than the homogeneous SFP. 
Moreover, to understand better the stabilization of the 
BO-modulated staggered flux wave-functions w.r.t the 
homogeneous one, we have plotted in Fig.HJthe difference 
in the respective kinetic energy, the exchange energy and 
the Coulomb energy of the SFPJ J and BO/J wave- 
functions. We conclude that the two wave- functions, al- 
though qualitatively similar (they both exhibit an un- 
derlying staggered flux pattern), bear quantitative dif- 
ferences: the staggered flux phase (slightly) better opti- 
mizes the kinetic energy whereas the bond-order wave- 
function (slightly) better optimizes the Coulomb and ex- 
change energies so that a small overall energy gain is in 
favor of the modulated phase. Therefore, we unambigu- 
ously conclude that, generically, bond-order modulations 
should spontaneously appear on top of the staggered flux 
pattern for moderate doping. 

Finally, we emphasize that the bond-order wave- 
function is not stabilized by the Coulomb repulsion alone 
(like for a usual electronic Wigner cristal) exhibiting co- 
existing bond order and (small) charge density wave. 
Moreover, the variational parameters in Eq. Q?JJ are 
found after minimizing the projected energy to be set to 
equal values on every site of the unit-cell. Let us also 
emphasize that the bond-order wave-function is not su- 
perconducting as proposed in some scenarios^. In the 
actual variational framework, we do not consider bond- 
order wave-function embedded in a sea of d-wave spin 
singlet pairs. 

In fact, we do not expect a bulk d-wave RVB state 
to be stable at large Coulomb repulsion (because of its 



very poor Coulomb energy) nor a bulk static checker- 
board SFP at too small Coulomb energy. However, for 
moderate Coulomb repulsion for which the d-wave RVB 
remains globally stable, sizeable regions of checkerboard 
SFP could be easily nucleated e.g. by defects. This is- 
sue will be addressed using renormalized MF theory in a 
future work. An extension of our VMC study with simul- 
taneous inhomogeneous bond-order and singlet pair order 
parameters (as required to treat such a problem) is diffi- 
cult and also left for a future work. Note also that low- 
energy dynamic fluctuations of checkerboard (and SFP) 
characters could also exist within the d-wave RVB state 
but this is beyond the scope of this paper. 

The properties of the BO/J staggered flux wave- 
function are summarized in Fig. [5] showing the real and 
imaginary parts of the measured hopping term (cfcj) 
between every nearest neighbor sites of our candidate 
BO j J wave- function. We also present the exchange 
term on each bonds of the lattice, and the local on- 
site charge density. We find that the bond-order wave- 
function has both (spin-spin) bond density wave and 
(small) charge density wave components. Nonetheless, 
the charge modulations are very small (the maximum 
charge deviation from the mean on-site charge is of the 
order of 2%) , and the charge density is a little bit 
larger in the center of the unit-cell. As expected, the 
SFPjJ has homogeneous hopping and exchange bonds 
within the unit-cell. Therefore, we conclude that af- 
ter projection the modulated variational wave-function 
differs quantitatively from the uniform one: the BO I J 
staggered flux wave-function is quite inhomogeneous (al- 
though with a very small charge modulation) leading to 
an increased magnetic energy gain while still preserving 
a competitive kinetic energy, a characteristic of the ho- 
mogeneous SFPj J wave-function. 



IV. CONCLUSION 

In conclusion, in this work we have investigated the 
t—J—V model using both mean-field calculations as well 
as more involved variational Monte-Carlo calculations. 
Both approaches have provided strong evidence that 
bond-order wave-functions (of underlying staggered flux 
character) are stabilized at zero and finite Coulomb re- 
pulsion for doping close to 1/8. In particular, variational 
Monte-Carlo calculations show that a bond modulation 
appears spontaneously on top of the staggered flux phase. 
This is in agreement with the work of Wang et al^S pre- 
dicting an instability of staggered flux type. We have 
also shown that the modulated and homogeneous SFP, 
although nearby in parameter space, are nevertheless sep- 
arated from each other by a small energy barrier. While 
both staggered flux wave- functions provide an optimal 
kinetic energy, the bond-modulated one exhibits a small 
extra gain of the exchange energy. On the other hand, 
a short range Coulomb repulsion favors both staggered 
flux wave-function w.r.t the d-wave RVB superconduc- 



9 




FIG. 9: (Color online) Local expectation values (a,b,c) of the kinetic and exchange energies of the projected BO j J wave- 
functions on each of the bonds within the unit-cell. Width of filled square symbols is proportional to the (a) real and (b) 
imaginary part of (cfcj), and (c) to the local exchange energy (Si • Sj). The sign of the imaginary part of the hopping bonds 
is according to the staggered flux pattern (arrows). The wave- function has small charge density variations (d), therefore we 
subtract the mean value n to the local density: size of circles are proportional to {m — n), and circles are open (filled) for 
negative (positive) sign. The biggest circle corresponds to an on-site charge deviation of 2%. All the above results are for Iq — 4 
and Vb = 5. 



tors and brings them close in energy. 

Finally, it would be interesting to study if the checker- 
board pattern could spontaneously appear in the vicinity 
of a vortex in the mixed phase of the cuprates. Such an 
issue could be addressed by studying the t— J— V model 
on a square lattice extending our variational scheme to 
include simultaneously nearest neighbor pairing and bond 
modulated staggered currents. It is expected that, while 
the pairing is suppressed in the vicinity of the vortex, the 
checkerboard pattern might be variationally stabilized in 



this region. 
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